****************************************************************************************** 
****************************************************************************************** 
****************************************************************************************** 
/*

THIS DO-FILE:

Gets variance-covariance matrix of the moments

*/
****************************************************************************************** 
****************************************************************************************** 
****************************************************************************************** 

clear all
set more off
set maxvar 10000

use "my_sample_long.dta", clear
merge 1:1 hhidpn wave using "imputed_kid_assets.dta", nogen
rename assets passets

gen 	agecohort=67+(wave-4)*2 if cohort==1
replace agecohort=73+(wave-4)*2 if cohort==2 // should be 72 but in biennial model the even numbers mess up
replace agecohort=77+(wave-4)*2 if cohort==3
replace agecohort=83+(wave-4)*2 if cohort==4 // should be 82...
replace agecohort=87+(wave-4)*2 if cohort==5

*turn age into a bi-yearly thing
gen agebi=2*ceil(age/2)

tab agecohort, gen(agecdum)

********************************
cap program drop momentboot
program prog_momentboot, eclass

	tempname ests
	mat `ests'=J(1,9+9+9+7+5+9+9+9+7+5+5+5+1+5,.)

	*kid assets
	forval a=1/9 {
		sum kassets if agecdum`a'==1 & cohort==1, det
		mat `ests'[1,`a']=r(p50)
	}
	forval a=1/9 {
		local b=`a'+3
		sum kassets if agecdum`b'==1 & cohort==2, det
		mat `ests'[1,9+`a']=r(p50)
	}
	forval a=1/9 {
		local b=`a'+5
		sum kassets if agecdum`b'==1 & cohort==3, det
		mat `ests'[1,18+`a']=r(p50)
	}
	forval a=1/7 {
		local b=`a'+8
		sum kassets if agecdum`b'==1 & cohort==4, det
		mat `ests'[1,27+`a']=r(p50)
	}	
	forval a=1/5 {
		local b=`a'+10
		sum kassets if agecdum`b'==1 & cohort==5, det
		mat `ests'[1,34+`a']=r(p50)
	}		
	*parent assets
	forval a=1/9 {
		sum passets if agecdum`a'==1 & cohort==1, det
		mat `ests'[1,39+`a']=r(p50)
	}
	forval a=1/9 {
		local b=`a'+3
		sum passets if agecdum`b'==1 & cohort==2, det
		mat `ests'[1,48+`a']=r(p50)
	}
	forval a=1/9 {
		local b=`a'+5
		sum passets if agecdum`b'==1 & cohort==3, det
		mat `ests'[1,57+`a']=r(p50)
	}
	forval a=1/7 {
		local b=`a'+8
		sum passets if agecdum`b'==1 & cohort==4, det
		mat `ests'[1,66+`a']=r(p50)
	}	
	forval a=1/5 {
		local b=`a'+10
		sum passets if agecdum`b'==1 & cohort==5, det
		mat `ests'[1,73+`a']=r(p50)
	}	
	
	*defining asset quintiles:
	*keeping it as all data cuz otherwise weight is too low... again the SD/SE question!
	cap drop apq
	xtile apq=passets, nq(5)

	*insurance
	forval x=1/5 {
		sum ltci if apq==`x'
		mat `ests'[1,78+`x']=r(mean)
	}

	*informal care
	forval x=1/5 {
		sum informaldummy if inlist(careneed,2,3) & apq==`x'
		mat `ests'[1,83+`x']=r(mean)
	}
	
	*medicaid
	qui sum medicaid
	mat `ests'[1,89]=r(mean)
	
	*kid full time
	forval x=1/5 {
		sum kidworkFT if cohort==`x'
		mat `ests'[1,89+`x']=r(mean)
	}	

	*mat list `ests'
	ereturn post `ests'
end
********************************

*run program 
prog_momentboot

*bootstrap
bootstrap _b, reps(1000) seed(712) cluster(hhidpn): prog_momentboot

mat vc_matrix=e(V)

mat list vc_matrix 
svmat vc_matrix

keep vc_matrix*

outsheet * using "variance_covariance_matrix.out", replace nonames
